home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Turnbull China Bikeride
/
Turnbull China Bikeride - Disc 1.iso
/
ARGONET
/
PD
/
MATHS
/
RLAB
/
RLAB125.ZIP
/
!RLaB
/
test
< prev
next >
Wrap
Text File
|
1995-11-18
|
73KB
|
3,194 lines
//
// Test Suite for RLaB.
// This test suite is supposed to test as much as possible, and still
// be runable on most platforms. This means we cannot do graphics or
// other platform specific stuff like piping. However, that is OK,
// since we are mostly interested in assuring the builder that RLaB
// will produce correct numerical answers.
// We use randsvd, which in turn uses RLaB's random number
// generator. We try to be carefull and seed the random number
// generator so that each user will get similar results (or at least
// similar inputs to the tests).
// Test the other style comments
1 + 2; // A simple statement with trailing comment.
# Optional RLaB comment style.
1 + 2; # A simple statement with trailing comment.
% Optional Matlab comment style.
1 + 2; % A simple statement with trailing comment.
srand (SEED = 10); // Seems to produce reasonable results
rand("default");
tic(); // Start timing the tests...
//
// Test Parameters and some functions we will need later
//
pi = 4.0*atan(1.0);
X = 3; // should be 3 (heuristic).
//
// Compute machine epsilon
//
epsilon = function()
{
eps = 1.0;
while((eps + 1.0) != 1.0)
{
eps = eps/2.0;
}
return 2*eps;
};
eps = epsilon();
eye = function( m , n )
{
if (!exist (n))
{
if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
new = zeros (m[1], m[2]);
N = min ([m[1], m[2]]);
else
if (class (m) == "string" || class (n) == "string") {
error ("eye(), string arguments not allowed");
}
if (max (size (m)) == 1 && max (size (n)) == 1)
{
new = zeros (m[1], n[1]);
N = min ([m[1], n[1]]);
else
error ("matrix arguments to eye() must be 1x1");
}
}
for(i in 1:N)
{
new[i;i] = 1.0;
}
return new;
};
symm = function( A )
{
return (A + A')./2;
};
//-------------------------------------------------------------------//
// Synopsis: Pascal matrix.
// Syntax: P = pascal ( N )
// Description:
// The Pascal matrix of order N: a symmetric positive definite
// matrix with integer entries taken from Pascal's triangle. The
// Pascal matrix is totally positive and its inverse has integer
// entries. Its eigenvalues occur in reciprocal pairs. COND(P)
// is approximately 16^N/(N*PI) for large N. PASCAL(N,1) is the
// lower triangular Cholesky factor (up to signs of columns) of
// the Pascal matrix. It is involutary (is its own
// inverse). PASCAL(N,2) is a transposed and permuted version of
// PASCAL(N,1) which is a cube root of the identity.
// References:
// R. Brawer and M. Pirovino, The linear algebra of the Pascal matrix,
// Linear Algebra and Appl., 174 (1992), pp. 13-23 (this paper
// gives a factorization of L = PASCAL(N,1) and a formula for the
// elements of L^k).
// S. Karlin, Total Positivity, Volume 1, Stanford University Press,
// 1968. (Page 137: shows i+j-1 choose j is TP (i,j=0,1,...).
// PASCAL(N) is a submatrix of this matrix.)
// M. Newman and J. Todd, The evaluation of matrix inversion programs,
// J. Soc. Indust. Appl. Math., 6(4):466--476, 1958.
// H. Rutishauser, On test matrices, Programmation en Mathematiques
// Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165,
// 1966, pp. 349-365. (Gives an integral formula for the
// elements of PASCAL(N).)
// J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
// Birkhauser, Basel, and Academic Press, New York, 1977, p. 172.
// H.W. Turnbull, The Theory of Determinants, Matrices, and Invariants,
// Blackie, London and Glasgow, 1929. (PASCAL(N,2) on page 332.)
// This file is a translation of pascal.m from version 2.0 of
// "The Test Matrix Toolbox for Matlab", described in Numerical
// Analysis Report No. 237, December 1993, by N. J. Higham.
//-------------------------------------------------------------------//
pascal = function ( n , k )
{
local(n, k)
if (!exist (k)) { k = 0; }
P = diag( (-1).^[0:n-1] );
P[; 1] = ones(n,1);
// Generate the Pascal Cholesky factor (up to signs).
for (j in 2:n-1)
{
for (i in j+1:n)
{
P[i;j] = P[i-1;j] - P[i-1;j-1];
}
}
if (k == 0)
{
P = P*P';
else if (k == 2) {
P = rot90(P,3);
if (n/2 == round(n/2)) { P = -P; }
}}
return P;
};
//-------------------------------------------------------------------//
// Synopsis: Random, orthogonal upper Hessenberg matrix.
// Syntax: H = ohess ( N )
// Description:
// H is an N-by-N real, random, orthogonal upper Hessenberg
// matrix. Alternatively, H = OHESS(X), where X is an arbitrary
// real N-vector (N > 1) constructs H non-randomly using the
// elements of X as parameters. In both cases H is constructed
// via a product of N-1 Givens rotations.
// Note: See Gragg (1986) for how to represent an N-by-N
// (complex) unitary Hessenberg matrix with positive subdiagonal
// elements in terms of 2N-1 real parameters (the Schur
// parametrization). This M-file handles the real case only and
// is intended simply as a convenient way to generate random or
// non-random orthogonal Hessenberg matrices.
// Reference:
// W.B. Gragg, The QR algorithm for unitary Hessenberg matrices,
// J. Comp. Appl. Math., 16 (1986), pp. 1-8.
// This file is a translation of ohess.m from version 2.0 of
// "The Test Matrix Toolbox for Matlab", described in Numerical
// Analysis Report No. 237, December 1993, by N. J. Higham.
//-------------------------------------------------------------------//
ohess = function ( x )
{
local (x)
global (pi)
if (any (imag (x))) { error("Parameter must be real."); }
n = max(size(x));
if (n == 1)
{
// Handle scalar x.
n = x;
x = rand(n-1, 1)*2*pi;
H = eye(n,n);
H[n;n] = sign(rand());
else
H = eye(n,n);
H[n;n] = sign(x[n]) + (x[n]==0); // Second term ensures H[n;n] nonzero.
}
for (i in n:2:-1)
{
// Apply Givens rotation through angle x(i-1).
theta = x[i-1];
c = cos(theta);
s = sin(theta);
H[ [i-1, i] ;] = [ c*H[i-1;]+s*H[i;] ;
-s*H[i-1;]+c*H[i;] ];
}
return H;
};
//-------------------------------------------------------------------//
// Synopsis: Random matrix with pre-assigned singular values.
// Syntax: R = randsvd (N, KAPPA, MODE, KL, KU)
// Description:
// R is a (banded) random matrix of order N with COND(A) = KAPPA
// and singular values from the distribution MODE.
// N may be a 2-vector, in which case the matrix is N(1)-by-N(2).
// Available types:
// MODE = 1: one large singular value,
// MODE = 2: one small singular value,
// MODE = 3: geometrically distributed singular values,
// MODE = 4: arithmetically distributed singular values,
// MODE = 5: random singular values with unif. dist. logarithm.
// If omitted, MODE defaults to 3, and KAPPA defaults to SQRT(1/EPS).
// If MODE < 0 then the effect is as for ABS(MODE) except that in the
// original matrix of singular values the order of the diagonal entries
// is reversed: small to large instead of large to small.
// KL and KU are the lower and upper bandwidths respectively; if they
// are omitted a full matrix is produced.
// If only KL is present, KU defaults to KL.
// Special case: if KAPPA < 0 then a random full symmetric positive
// definite matrix is produced with COND(A) = -KAPPA and
// eigenvalues distributed according to MODE.
// KL and KU, if present, are ignored.
// This routine is similar to the more comprehensive Fortran routine xLATMS
// in the following reference:
// J.W. Demmel and A. McKenney, A test matrix generation suite,
// LAPACK Working Note #9, Courant Institute of Mathematical Sciences,
// New York, 1989.
// This file is a translation of randsvd.m from version 2.0 of
// "The Test Matrix Toolbox for Matlab", described in Numerical
// Analysis Report No. 237, December 1993, by N. J. Higham.
// Dependencies
# rfile bandred
# rfile qmult
//-------------------------------------------------------------------//
randsvd = function (n, kappa, mode, kl, ku)
{
local (n, kappa, mode, kl, ku)
if (!exist (kappa)) { kappa = sqrt(1/eps); }
if (!exist (mode)) { mode = 3; }
if (!exist (kl)) { kl = n-1; } // Full matrix.
if (!exist (ku)) { ku = kl; } // Same upper and lower bandwidths.
if (abs(kappa) < 1) {
error("Condition number must be at least 1!");
}
posdef = 0;
if (kappa < 0) // Special case.
{
posdef = 1;
kappa = -kappa;
}
p = min(n);
m = n[1]; // Parameter n specifies dimension: m-by-n.
n = n[max(size(n))];
if (p == 1) // Handle case where A is a vector.
{
rand("normal", -10, 10);
A = rand(m, n);
A = A/norm(A, "2");
return A;
}
j = abs(mode);
// Set up vector sigma of singular values.
if (j == 3)
{
factor = kappa^(-1/(p-1));
sigma = factor.^[0:p-1];
else if (j == 4) {
sigma = ones(p,1) - (0:p-1)'/(p-1)*(1-1/kappa);
else if (j == 5) { // In this case cond(A) <= kappa.
rand("uniform", 0, 1)
sigma = exp( -rand(p,1)*log(kappa) );
else if (j == 2) {
sigma = ones(p,1);
sigma[p] = 1/kappa;
else if (j == 1) {
sigma = ones(p,1)./kappa;
sigma[1] = 1;
}}}}}
// Convert to diagonal matrix of singular values.
if (mode < 0) {
sigma = sigma[p:1:-1];
}
sigma = diag(sigma);
if (posdef) // Handle special case.
{
Q = qmult(p);
A = Q'*sigma*Q;
A = (A + A')/2; // Ensure matrix is symmetric.
return A;
}
if (m != n)
{
sigma[m; n] = 0; // Expand to m-by-n diagonal matrix.
}
if (kl == 0 && ku == 0) // Diagonal matrix requested - nothing more to do.
{
A = sigma;
return A;
}
// A = U*sigma*V, where U, V are random orthogonal matrices from the
// Haar distribution.
A = qmult(sigma');
A = qmult(A');
if (kl < n-1 || ku < n-1) // Bandwidth reduction.
{
A = bandred(A, kl, ku);
}
rand("default");
return A;
};
//-------------------------------------------------------------------//
// Synopsis: Band reduction by two-sided unitary transformations.
// Syntax: bandred ( A , KL, KU )
// Description:
// bandred(A, KL, KU) is a matrix unitarily equivalent to A with
// lower bandwidth KL and upper bandwidth KU (i.e. B(i,j) = 0 if
// i > j+KL or j > i+KU). The reduction is performed using
// Householder transformations. If KU is omitted it defaults to
// KL.
// Called by randsvd.
// This is a `standard' reduction. Cf. reduction to bidiagonal
// form prior to computing the SVD. This code is a little
// wasteful in that it computes certain elements which are
// immediately set to zero!
//
// Reference:
// G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
// Johns Hopkins University Press, Baltimore, Maryland, 1989.
// Section 5.4.3.
// This file is a translation of bandred.m from version 2.0 of
// "The Test Matrix Toolbox for Matlab", described in Numerical
// Analysis Report No. 237, December 1993, by N. J. Higham.
// Dependencies
# rfile house
//-------------------------------------------------------------------//
bandred = function ( A , kl , ku )
{
local (A, kl, ku)
if (!exist (ku)) { ku = kl; else ku = ku; }
if (kl == 0 && ku == 0) {
error ("You''ve asked for a diagonal matrix. In that case use the SVD!");
}
// Check for special case where order of left/right transformations matters.
// Easiest approach is to work on the transpose, flipping back at the end.
flip = 0;
if (ku == 0)
{
A = A';
temp = kl; kl = ku; ku = temp; flip = 1;
}
m = A.nr; n = A.nc;
for (j in 1 : min( min(m, n), max(m-kl-1, n-ku-1) ))
{
if (j+kl+1 <= m)
{
ltmp = house(A[j+kl:m;j]);
beta = ltmp.beta; v = ltmp.v;
temp = A[j+kl:m;j:n];
A[j+kl:m;j:n] = temp - beta*v*(v'*temp);
A[j+kl+1:m;j] = zeros(m-j-kl,1);
}
if (j+ku+1 <= n)
{
ltmp = house(A[j;j+ku:n]');
beta = ltmp.beta; v = ltmp.v;
temp = A[j:m;j+ku:n];
A[j:m;j+ku:n] = temp - beta*(temp*v)*v';
A[j;j+ku+1:n] = zeros(1,n-j-ku);
}
}
if (flip) {
A = A';
}
return A;
};
//-------------------------------------------------------------------//
// Synopsis: Lehmer matrix - symmetric positive definite.
// Syntax: A = lehmer ( N )
// Description:
// A is the symmetric positive definite N-by-N matrix with
// A(i,j) = i/j for j >= i.
// A is totally nonnegative. INV(A) is tridiagonal, and explicit
// formulas are known for its entries.
// N <= COND(A) <= 4*N*N.
// References:
// M. Newman and J. Todd, The evaluation of matrix inversion
// programs, J. Soc. Indust. Appl. Math., 6 (1958), pp. 466-476.
// Solutions to problem E710 (proposed by D.H. Lehmer): The inverse
// of a matrix, Amer. Math. Monthly, 53 (1946), pp. 534-535.
// J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
// Birkhauser, Basel, and Academic Press, New York, 1977, p. 154.
// This file is a translation of lehmer.m from version 2.0 of
// "The Test Matrix Toolbox for Matlab", described in Numerical
// Analysis Report No. 237, December 1993, by N. J. Higham.
//-------------------------------------------------------------------//
lehmer = function ( n )
{
local (n)
global (tril)
A = ones(n,1)*(1:n);
A = A./A';
A = tril(A) + tril(A,-1)';
return A;
};
//-------------------------------------------------------------------//
// Synopsis: Pre-multiply by random orthogonal matrix.
// Syntax: B = qmult ( A )
// Description:
// B is Q*A where Q is a random real orthogonal matrix from the
// Haar distribution, of dimension the number of rows in
// A. Special case: if A is a scalar then QMULT(A) is the same as
// qmult(eye(a))
// Called by RANDSVD.
// Reference:
// G.W. Stewart, The efficient generation of random
// orthogonal matrices with an application to condition estimators,
// SIAM J. Numer. Anal., 17 (1980), 403-409.
// This file is a translation of qmult.m from version 2.0 of
// "The Test Matrix Toolbox for Matlab", described in Numerical
// Analysis Report No. 237, December 1993, by N. J. Higham.
//-------------------------------------------------------------------//
qmult = function ( A )
{
local (A)
n = A.nr; m = A.nc;
// Handle scalar A.
if (max(n,m) == 1)
{
n = A;
A = eye(n,n);
}
d = zeros(n,n);
for (k in n-1:1:-1)
{
// Generate random Householder transformation.
rand("normal", 0, 1);
x = rand(n-k+1,1);
s = norm(x, "2");
sgn = sign(x[1]) + (x[1]==0); // Modification for sign(1)=1.
s = sgn*s;
d[k] = -sgn;
x[1] = x[1] + s;
beta = s*x[1];
// Apply the transformation to A.
y = x'*A[k:n;];
A[k:n;] = A[k:n;] - x*(y/beta);
}
// Tidy up signs.
for (i in 1:n-1)
{
A[i;] = d[i]*A[i;];
}
A[n;] = A[n;]*sign(rand());
B = A;
rand("default");
return B;
};
//-------------------------------------------------------------------//
// Synopsis: Create a Householder matrix.
// Syntax: house ( X )
// If HOUSE(x), which returns a list containing elements `v' and
// `beta', then H = EYE - beta*v*v' is a Householder matrix such
// that Hx = -sign(x(1))*norm(x)*e_1.
// NB: If x = 0 then v = 0, beta = 1 is returned.
// x can be real or complex.
// sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0).
// Theory: (textbook references Golub & Van Loan 1989, 38-43;
// Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50).
// Hx = y: (I - beta*v*v')x = -s*e_1.
// Must have |s| = norm(x), v = x+s*e_1, and
// x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
// So take s = sign(x(1))*norm(x) (which avoids cancellation).
// v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
// = 2*norm(x)*(norm(x) + |x(1)|).
// References:
// G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
// Johns Hopkins University Press, Baltimore, Maryland, 1989.
// G.W. Stewart, Introduction to Matrix Computations, Academic Press,
// New York, 1973,
// J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
// Press, 1965.
// This file is a translation of house.m from version 2.0 of
// "The Test Matrix Toolbox for Matlab", described in Numerical
// Analysis Report No. 237, December 1993, by N. J. Higham.
//-------------------------------------------------------------------//
house = function ( x )
{
local (x)
m = x.nr; n = x.nc;
if (n > 1) { error ("Argument must be a column vector."); }
s = norm(x,"2") * (sign(x[1]) + (x[1]==0)); // Modification for sign(0)=1.
v = x;
if (s == 0) // Quit if x is the zero vector.
{
beta = 1;
return << beta = beta ; v = v >>;
}
v[1] = v[1] + s;
beta = 1/(s'*v[1]); // NB the conjugated s.
// beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real.
// But beta as above can be non-real (due to rounding) only
// when x is complex.
return << beta = beta ; v = v >>;
};
//-------------------------------------------------------------------//
// Syntax: tril ( A )
// tril ( A , K )
// Description:
// tril(x) returns the lower triangular part of A.
// tril(A,K) returns the elements on and below the K-th diagonal of
// A.
// K = 0: main diagonal
// K > 0: above the main diag.
// K < 0: below the main diag.
// See Also: triu
//-------------------------------------------------------------------//
tril = function(x, k)
{
local(x, k)
if (!exist (k)) { k = 0; }
nr = x.nr; nc = x.nc;
if(k > 0)
{
if (k > (nc - 1)) { error ("tril: invalid value for k"); }
else
if (abs (k) > (nr - 1)) { error ("tril: invalid value for k"); }
}
y = zeros(nr, nc);
for(i in max( [1,1-k] ):nr) {
j = 1:min( [nc, i+k] );
y[i;j] = x[i;j];
}
return y;
};
//-------------------------------------------------------------------//
// Syntax: triu ( A )
// triu ( A , K )
// Description:
// triu(x) returns the upper triangular part of A.
// tril(x; k) returns the elements on and above the k-th diagonal of
// A.
// K = 0: main diagonal
// K > 0: above the main diag.
// K < 0: below the main diag.
// See Also: tril
//-------------------------------------------------------------------//
triu = function(x, k)
{
local(x, k)
if (!exist (k)) { k = 0; }
nr = x.nr; nc = x.nc;
if(k > 0)
{
if (k > (nc - 1)) { error ("triu: invalid value for k"); }
else
if (abs (k) > (nr - 1)) { error ("triu: invalid value for k"); }
}
y = zeros(nr, nc);
for(j in max( [1,1+k] ):nc) {
i = 1:min( [nr, j-k] );
y[i;j] = x[i;j];
}
return y;
};
//
//-------------------- Test Relational Expressions -------------------
//
printf("\tstart scalar tests...\n");
printf("\tstart relational tests...\n");
// SCALAR CONSTANTS (REAL)
if( !(1<2) ) { error(); }
if( !(1<=2) ) { error(); }
if( 1>2 ) { error(); }
if( 1>=2 ) { error(); }
if( 1==2 ) { error(); }
if( !(1!=2) ) { error(); }
if( !1 ) { error(); }
if( !!!1) { error(); }
if( !([1]<[2]) ) { error(); }
if( !([1]<=[2]) ) { error(); }
if( [1]>[2] ) { error(); }
if( [1]>=[2] ) { error(); }
if( [1]==[2] ) { error(); }
if( !([1]!=[2]) ) { error(); }
if( ![1] ) { error(); }
if( !!![1]) { error(); }
// SCALAR CONSTANTS (COMPLEX)
if( !(1+2i<2+3i) ) { error(); }
if( !(1+2i<=2+3i) ) { error(); }
if( 1+2i>2+3i ) { error(); }
if( 1+2i>=2+3i ) { error(); }
if( 1+2i==2+3i ) { error(); }
if( !(1+2i!=2+3i) ) { error(); }
if( !1+2i ) { error(); }
if( !!!1+2i) { error(); }
if( !([1+2i]<[2+3i]) ) { error(); }
if( !([1+2i]<=[2+3i]) ) { error(); }
if( [1+2i]>[2+3i] ) { error(); }
if( [1+2i]>=[2+3i] ) { error(); }
if( [1+2i]==[2+3i] ) { error(); }
if( !([1+2i]!=[2+3i]) ) { error(); }
if( ![1+2i] ) { error(); }
if( !!![1+2i]) { error(); }
// SCALAR ENTITIES (REAL)
a=1;b=2;
if( !(a<b) ) { error(); }
if( !(a<=b) ) { error(); }
if( a>b ) { error(); }
if( a>=b ) { error(); }
if( a==b ) { error(); }
if( !(a!=b) ) { error(); }
if( !a ) { error(); }
if( !!!a) { error(); }
if( !([a]<[b]) ) { error(); }
if( !([a]<=[b]) ) { error(); }
if( [a]>[b] ) { error(); }
if( [a]>=[b] ) { error(); }
if( [a]==[b] ) { error(); }
if( !([a]!=[b]) ) { error(); }
if( ![a] ) { error(); }
if( !!![a]) { error(); }
// SCALAR ENTITIES (COMPLEX)
a=1+2i;b=2+3i;
if( !(a<b) ) { error(); }
if( !(a<=b) ) { error(); }
if( a>b ) { error(); }
if( a>=b ) { error(); }
if( a==b ) { error(); }
if( !(a!=b) ) { error(); }
if( !a ) { error(); }
if( !!!a) { error(); }
if( !([a]<[b]) ) { error(); }
if( !([a]<=[b]) ) { error(); }
if( [a]>[b] ) { error(); }
if( [a]>=[b] ) { error(); }
if( [a]==[b] ) { error(); }
if( !([a]!=[b]) ) { error(); }
if( ![a] ) { error(); }
if( !!![a]) { error(); }
x = rand(4,4);
if (! all (all (-2 < x))) { "-2 < x" error (); }
if (! all (all (-2 <= x))) { "-2 <= x" error (); }
if (! all (all ( 2 > x))) { "2 > x" error (); }
if (! all (all ( 2 >= x))) { "2 >= x" error (); }
if (! all (all (rand(4,4) > -2))) { error (); }
if (! all (all (rand(4,4) >= -2))) { error (); }
if (! all (all (rand(4,4) < 2))) { error (); }
if (! all (all (rand(4,4) <= 2))) { error (); }
if (! all (all (rand (4,4) > -rand (4,4)))) { error (); }
if (! all (all (rand (4,4) >= -rand (4,4)))) { error (); }
if (! all (all (-rand (4,4) < rand (4,4)))) { error (); }
if (! all (all (-rand (4,4) <= rand (4,4)))) { error (); }
//------------------------- Test REAL SCALARS ------------------------
//
// CONSTANTS
// Addition
if(1+2 != 3) { error(); }
// Subtraction
if(1-2 != -1) { error(); }
// Multiply
if(1*2 != 2) { error(); }
// Divide
if(1/2 != 0.5) { error(); }
// Power
if(2^2 != 4) { error(); }
if(4^0 != 1) { error(); }
// Unary Minus
if(2-3 != -1) { error(); }
//
// ENTITIES
//
a = 1; b = 2; c = 3; d = 0.5;
// Addition
if(a+b != c) { error(); }
// Subtraction
if(a-b != -a) { error(); }
// Multiply
if(a*b != b) { error(); }
// Divide
if(a/b != d) { error(); }
// Power
if(b^b != b*b) { error(); }
if(b^0 != 1) { error (); }
// Unary Minus
if(b-c != -a) { error(); }
//
// ENTITIES & CONSTANTS
//
if(a+2 != c) { error(); }
// Subtraction
if(1-b != -a) { error(); }
// Multiply
if(a*2 != 2) { error(); }
// Divide
if(1/b != d) { error(); }
// Power
if(2^b != b*b) { error(); }
// Unary Minus
if(b-3 != -a) { error(); }
//
//------------------------Test COMPLEX SCALARS -------------------------
//
// CONSTANTS
if(sqrt(-1) != 1i) { error(); }
// Addition
if((1+2i)+(2+3i) != (3+5i)) { error(); }
// Subtraction
if((1+2i)-(3+4i) != (-2-2i)) { error(); }
// Multiply
if((1+2i)*(3+4i) != (-5+10i)) { error(); }
// Divide
if((1+2i)/(3-4i) != (-.2+.4i)) { error(); }
// Power
// Precision problems prevent us from testing these. Have to
// be checked by hand.
// (1+2i)^2 = -3 + 4i
// (1+2i)^.5 = 1.272 + 7.862e-1i
// if((1+2i)^2 != (-3+4i)) { error(); }
// if((1+2i)^10 != (237+3116i)) { error(); }
// Unary Minus
if(-(1+2i) != -1-2i) { error(); }
//
// ENTITIES
//
a = 1+2i; b = 3+4i; c = 4+6i;
if(a+b != c) { error(); }
// Subtraction
if(a-b != -2-2i) { error(); }
// Multiply
if(a*b != -5+10i) { error(); }
// Divide
//if(a/(3-4i) != -.2+.4i) { error(); }
// Power
// if(b^b != b*b) { error(); }
// Unary Minus
if(-a != -1-2i) { error(); }
//
// ENTITIES & CONSTANTS
//
if(a+(3+4i) != c) { error(); }
// Subtraction
if((1+2i)-b != -2-2i) { error(); }
// Multiply
if(a*(3+4i) != -5+10i) { error(); }
// Divide
//if(a/(3-4i) != -.2+.4i) { error(); }
// Power
//if(b^b != b*b) { error(); }
// Unary Minus
if(-(1+2i) != -1-2i) { error(); }
//
// String - Numerical Equalities
//
if ((1 == "1")) { error(); }
if (([1] == "1")) { error(); }
if (("1" == 1)) { error(); }
if (("1" == [1])) { error(); }
if (rand(3,3) == "str") { error(); }
if ("str" == rand(3,3)) { error(); }
if (!any (any ((["1", "2"; "3", "4"] == "1") == [1,0;0,0]))) { error(); }
//
//----------------------- Test REAL MATRICES ---------------------------
//
printf("\tstart matrix tests...\n");
printf("\t\treal-matrices\n");
// Read in test matrices
//
read("test_input");
//
// Matrix construction
//
if(any([1;2;3] != [1,2,3]')) {
error();
}
if(any (any (m0 != zeros(2,2)))) { error(); }
if(any (any (m1 != 1+zeros(2,2)))) { error(); }
if(any (any (m2 != [1,2;3,4]))) { error(); }
if(any (any (m3 != [1+2i,2+3i;3+4i,5+6i]))) { error(); }
if(any (any ([1,2;3+0i,4+0i] != m2))) { error(); }
if(any (any ([m2] != m2))) { error(); }
//
// Matrix indexing
//
p = pascal(6);
if (!all (all (p[ [1:3] ; ] == p[ [1:3]' ; ]))) { error (); }
if (!all (all (p[ ; [1:3] ] == p[ ; [1:3]' ]))) { error (); }
if (!all (all (p[ [2:6] ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
if (!all (all (p[ [2:6]' ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
if (!all (all (p[ [6:12]' ] == p[ [6:12] ]))) { error (); }
//
// Sub-Matrix promotion
//
if(m2[2;2] != 4) { error(); }
if(any(m2[2;] != [3,4])) { error(); }
if(any(m2[;2] != [2,4]')) { error(); }
i=2;j=1;
if(m2[i;j] != 3) { error(); }
i=1;j=2;
if(m2[i;j] != 2) { error(); }
m = [1,2,3;4,5,6;7,8,9];
if(any(m[1;1,2] != [1,2]))
{
error();
}
if(any(m[1,2;1] != [1;4]))
{
error();
}
if(any (any (m[1,2;1,2] != [1,2;4,5])))
{
error();
}
//
if(m3[2;2] != (5+6i)) { error(); }
if(any(m3[2;] != [3+4i,5+6i])) { error(); }
if(any(m3[;2] != conj([2+3i,5+6i]'))) { error(); }
//
// Automatic creation, extension
//
if(any (any ((mm[3;3]=10) != [0,0,0;0,0,0;0,0,10]))) { error(); }
a=[1,2,3;4,5,6];
if(any (any ((a[3;1]=10) != [1,2,3;4,5,6;10,0,0]))) { error(); }
a=[1,2;3,4];
if(any (any ((a[3,4;3,4]=[5,6;7,8]) != [1,2,0,0;3,4,0,0;0,0,5,6;0,0,7,8])))
{
error();
}
mmm[2;] = a[1;];
//
// Matrix binary operations
//
a = m2; b = [5,6;7,8];
if(any (any (a+a != [2,4;6,8]))) { error(); }
if(any (any (a-a != zeros(2,2)))) { error(); }
if(any (any (2+a != [3,4;5,6]))) { error(); }
if(any (any (2-a != [1,0;-1,-2]))) { error(); }
if(any (any (a-2 != [-1,0;1,2]))) { error(); }
if(any (any (2*a != [2,4;6,8]))) { error(); }
if(any (any ((a./2 != [0.5,1;1.5,2])))) { error(); }
if(any (any (a*a != [7,10;15,22]))) { error(); }
if(any (any (a*a*a != [37,54;81,118]))) { error(); }
if(any (any (a .* a != [1,4;9,16]))) { error(); }
if(any (any (a./a != [1,1;1,1]))) { error(); }
if(any (any (a' != [1,3;2,4]))) { error(); }
if(any(any(rand(3,3)^0 != eye(3,3)))) { error(); }
if(any(any(rand(3,3).^0 != ones(3,3)))) { error(); }
if(any(any(rand(1,3).^0 != ones(1,3)))) { error(); }
if(any(any(rand(3,1).^0 != ones(3,1)))) { error(); }
if(any(any(1.^zeros(3,1) != ones(3,1)))) { error(); }
if(any(any(1.^zeros(1,3) != ones(1,3)))) { error(); }
if(any ([1;2;3]+[4;5;6] != [5;7;9]))
{
error();
}
if(any ([1;2;3]-[4;5;6] != [-3;-3;-3]))
{
error();
}
if(any ([2;2;2] ./ [1;1;1] != [2;2;2]))
{
error();
}
if(any ([1;2;3] .* [4;5;6] != [4;10;18]))
{
error();
}
if (type (1^.33333) != "real") { error (); }
if (type (1^(1/3)) != "real") { error (); }
if (type ([1]^.33333) != "real") { error (); }
if (type (1^[.33333]) != "real") { error (); }
if (type ([1]^[.33333]) != "real") { error (); }
//
// Test row-wise matrix addition
//
a = [1,2,3]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
if (!all (all (a .+ b == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
if (!all (all (b .+ a == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
a = [1,2,3] + [1,2,3]*1i;
b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
c = [2,4,6;5,7,9;8,10,12;11,13,15] + [2,4,6;5,7,9;8,10,12;11,13,15]*1i;
if (!all (all (a .+ b == c))) { error (); }
if (!all (all (b .+ a == c))) { error (); }
printf("\t\tpassed matrix row-wise add test...\n");
//
// Test row-wise matrix subtraction
//
a = [1,1,1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
if (!all (all (b .- a == [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
a = [1,1,1] + [1,1,1]*1i;
b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
if (!all (all (a .- b == -c))) { error (); }
if (!all (all (b .- a == c))) { error (); }
printf("\t\tpassed matrix row-wise subtraction test...\n");
//
// Test col-wise matrix addition
//
a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
if (!all (all (a .+ b == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
if (!all (all (b .+ a == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
a = [1;1;1;1] + [1;1;1;1]*1i;
b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
c = [2,3,4;5,6,7;8,9,10;11,12,13] + [2,3,4;5,6,7;8,9,10;11,12,13]*1i;
if (!all (all (a .+ b == c))) { error (); }
if (!all (all (b .+ a == c))) { error (); }
printf("\t\tpassed matrix col-wise add test...\n");
//
// Test col-wise matrix subtraction
//
a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
if (!all (all (b .- a == [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
a = [1;1;1;1] + [1;1;1;1]*1i;
b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
if (!all (all (a .- b == -c))) { error (); }
if (!all (all (b .- a == c))) { error (); }
printf("\t\tpassed matrix col-wise subtraction test...\n");
a = [1,2,3];
b = [1,2,3;4,5,6;7,8,9];
c = [1,4,9;4,10,18;7,16,27];
if (!all (all (a .* b == c))) { error (); }
if (!all (all (b .* a == c))) { error (); }
za = a + rand (size (a))*1j;
zb = b + rand (size (b))*1j;
if (!all (all (za .* zb == [za;za;za] .* zb))) { error (); }
if (!all (all (zb .* za == zb .* [za;za;za]))) { error (); }
if (!all (all (a .* zb == [a;a;a] .* zb))) { error (); }
if (!all (all (zb .* a == zb .* [a;a;a]))) { error (); }
if (!all (all (za .* b == [za;za;za] .* b))) { error (); }
if (!all (all (b .* za == b .* [za;za;za]))) { error (); }
printf("\t\tpassed matrix row-wise multiplication test...\n");
a = [1,2,3];
b = [1,2,3;4,6,6;7,8,9];
c = [1,1,1;4,3,2;7,4,3];
if (!all (all (b ./ a == c))) { error (); }
if (!all (all ([a;a;a] ./ b == a ./ b))) { error (); }
if (!all (all (b ./ [a;a;a] == b ./ a))) { error (); }
za = a + rand (size (a))*1j;
zb = b + rand (size (b))*1j;
if (!all (all ([za;za;za] ./ zb == za ./ zb))) { error (); }
if (!all (all (zb ./ [za;za;za] == zb ./ za))) { error (); }
if (!all (all ([a;a;a] ./ zb == a ./ zb))) { error (); }
if (!all (all (zb ./ [a;a;a] == zb ./ a))) { error (); }
if (!all (all ([za;za;za] ./ b == za ./ b))) { error (); }
if (!all (all (b ./ [za;za;za] == b ./ za))) { error (); }
printf("\t\tpassed matrix row-wise division test...\n");
a = [1;2;3];
b = [1,2,3;4,5,6;7,8,9];
if (!all (all (a .* b == [a,a,a] .* b))) { error (); }
if (!all (all (b .* a == b .* [a,a,a]))) { error (); }
za = a + rand (size (a))*1j;
zb = b + rand (size (b))*1j;
if (!all (all (za .* zb == [za,za,za] .* zb))) { error (); }
if (!all (all (zb .* za == zb .* [za,za,za]))) { error (); }
if (!all (all (za .* b == [za,za,za] .* b))) { error (); }
if (!all (all (b .* za == b .* [za,za,za]))) { error (); }
if (!all (all (a .* zb == [a,a,a] .* zb))) { error (); }
if (!all (all (zb .* a == zb .* [a,a,a]))) { error (); }
printf("\t\tpassed matrix column-wise multiplication test...\n");
a = [1;2;3];
b = [1,2,3;4,6,6;7,8,9];
if (!all (all ([a,a,a] ./ b == a ./ b))) { error (); }
if (!all (all (b ./ [a,a,a] == b ./ a))) { error (); }
za = a + rand (size(a))*1j;
zb = b + rand (size(b))*1j;
if (!all (all ([za,za,za] ./ zb == za ./ zb))) { error (); }
if (!all (all (zb ./ [za,za,za] == zb ./ za))) { error (); }
if (!all (all ([za,za,za] ./ b == za ./ b))) { error (); }
if (!all (all (b ./ [za,za,za] == b ./ za))) { error (); }
if (!all (all ([a,a,a] ./ zb == a ./ zb))) { error (); }
if (!all (all (zb ./ [a,a,a] == zb ./ a))) { error (); }
printf("\t\tpassed matrix column-wise division test...\n");
//
//--------------------- Test COMPLEX MATRICES --------------------------
//
// Automatic creation, extension
//
printf("\t\tcomplex-matrices\n");
if(any (any ((mm[3;3]=10+10i) != [0,0,0;0,0,0;0,0,10+10i]))) { error(); }
a=[1,2,3;4,5,6];
if(any (any ((a[3;1]=10+10i) != [1,2,3;4,5,6;10+10i,0,0]))) { error(); }
//
a = m3;
if(any (any (a+a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
if(any (any (a-a != zeros(2,2)))) { error(); }
if(any (any (2+a != [3+2i,4+3i;5+4i,7+6i]))) { "2+a" error(); }
if(any (any (2-a != [1-2i,0-3i;-1-4i,-3-6i]))) { error(); }
if(any (any (a-2 != [-1+2i,0+3i;1+4i,3+6i]))) { error(); }
if(any (any (2*a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
if(any (any (a./2 != [.5+1i,1+1.5i;1.5+2i,2.5+3i]))) { error(); }
if(any (any (a*a != [-9+21i,-12+34i;-14+48i,-17+77i]))) { error(); }
if(any (any (a*a*a != [-223+57i,-345+113i;-469+183i,-719+337i]))) { error(); }
if(any (any (a .* a != [-3+4i,-5+12i;-7+24i,-11+60i]))) { error(); }
//
// The following test may not work on some machines
//
if(any (any (a./a != [1,1;1,1]))) {
printf("\t\t***complex division inaccuracy, check manually***\n");
}
if(any (any (a' != conj([1+2i,3+4i;2+3i,5+6i])))) { error(); }
//
//--------------------- Test NULL MATRICES -------------------------
//
printf("\t\tnull-matrices\n");
// Create a NULL matrix
mnull = [];
// Test it with SCALARS
if( any([1,mnull] != 1)) {
error();
}
if( any([mnull,1] != 1)) {
error();
}
// Test with MATRIX construction
m = [1,2;3,4;5,6];
if( any([mnull;1] != [1])) {
error();
}
if( any([1;mnull] != [1])) {
error();
}
if( any([mnull;1,2,3] != [1,2,3])) {
error();
}
if( any([1,2,3;mnull] != [1,2,3])) {
error();
}
if(any (any ([mnull,m] != m))) {
error();
}
if(any (any ([m,mnull] != m))) {
error();
}
if(any (any ([mnull;m] != m))) {
error();
}
if(any (any ([m;mnull] != m))) {
error();
mnull = matrix();
mnull[1] = [1];
}
//--------------------- Test Matrix Multiply --------------------------
i = sqrt(-1);
a = [1,2,3;4,5,6;7,8,9];
b = [4,5,6;7,8,9;10,11,12];
c = [ 48, 54, 60 ;
111, 126, 141 ;
174, 198, 222 ];
if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
az = a + b*i;
bz = b + a*i;
cz = [-18+141*i , -27+162*i , -36+183*i ;
9+240*i , 0+279*i , -9+318*i ;
36+339*i , 27+396*i , 18+453*i ];
czz = [ 48+30*i , 54+36*i , 60+42*i ;
111+66*i , 126+81*i , 141+96*i ;
174+102*i, 198+126*i , 222+150*i ];
czzz = [ 48+111*i , 54+126*i , 60+141*i ;
111+174*i , 126+198*i , 141+222*i ;
174+237*i , 198+270*i , 222+303*i ];
if (any (any (cz != az*bz))) { error ("failed Complex-Complex Multiply"); }
if (any (any (czz != a*bz))) { error ("failed Real-Complex Multiply"); }
if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
a = [a,a];
b = [b;b];
c = [ 96 , 108 , 120 ;
222 , 252 , 282 ;
348 , 396 , 444 ];
if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
az = [az,az];
bz = [bz;bz];
cz = [ -36+282*i , -54+324*i , -72+366*i ;
18+480*i , 0+558*i , -18+636*i ;
72+678*i , 54+792*i , 36+906*i ];
czz = [ 96+60*i , 108+72*i , 120+84*i ;
222+132*i , 252+162*i , 282+192*i ;
348+204*i , 396+252*i , 444+300*i ];
czzz = [ 96+222*i , 108+252*i , 120+282*i ;
222+348*i , 252+396*i , 282+444*i ;
348+474*i , 396+540*i , 444+606*i ];
if (any (any (cz != az*bz))) { error ("failed Complex-Complex Multiply"); }
if (any (any (czz != a*bz))) { error ("failed Real-Complex Multiply"); }
if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
printf("\t\tpassed matrix multiply test...\n");
//--------------------- Test STRING MATRICES --------------------------
//
printf("\t\tstring-matrices\n");
sm = ["s1","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"];
if(sm[1] != "s1") { error(); }
if( sm[1;3] != "sm3" ) { error(); }
if(any(sm[2,3;3] != ["xxx3";"yyy3"]) ) { error(); }
if(any (any ((sm[1;1]="xx")!=["xx","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"])))
{
error();
}
if( ((strm[1] = "strm")[1]) != "strm" ) { error(); }
// Test string-matrix add functionality
sm = sm[1,2;1,2];
if (any (any (("1_"+sm+"_2") != ["1_xx_2","1_sm2_2";"1_x1_2","1_x2_2"]))) {error();}
if ("c"+["1"] != "c1") { error (); }
if (["c"]+"1" != "c1") { error (); }
// An append test...
bs = "b";
bm = "b";
xsm = ["a", bm, bs ];
if (!all(xsm == ["a", "b", "b"])) { error (); }
printf("\tpassed matrix test...\n");
//
//---------------------------- List Tests --------------------------
//
// List creation
listest = << << 11; 12 >>; << 21; 22>> >>;
if( listest.[1].[2] != 12 ) { error(); }
if(any(<<a=10;b=1:4;c=[1,2,3;4,5,6;7,8,9]>>.b != [1,2,3,4])) { error(); }
mlist.[0] = m;
if(any(any(mlist.[0] != m))) { error(); }
// Test list functions...
listtest.fun = function ( a ) { return sin(a); };
if (!(listtest.fun(0.5) == sin(0.5))) { error() ; }
// Test open-list assignment...
</v;d/> = eig(rand(3,3));
printf("\tpassed list test...\n");
//
// Reset random number generator seed...
//
rand("default");
//
//-------------------------- Test printf () --------------------------
//
sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", 5,3,2, 8,7,3, "string", 3, 4, 1234e-2);
if (!(tmp == " 002 0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", [5],[3],[2], [8],[7],[3], ...
["string"], [3], [4], [1234e-2]);
if (!(tmp == " 002 0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", [5,1],[3,2],[2,2], [8],[7],[3], ...
["string"], [3], [4], [1234e-2,4]);
if (!(tmp == " 002 0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", [5+2i,1],[3,2+4i],[2,2], [8],[7],[3], ...
["string"], [3+2i], [4], [1234e-2+12j,4]);
if (!(tmp == " 002 0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
printf("\tpassed sprintf test...\n");
//
//------------------------- Test strtod() ----------------------------
//
if (123.456 != strtod ("123.456")) { error (); }
if (!all (all ([1,2;3,4] == strtod (["1","2";"3","4"]))))
{ error (); }
printf("\tpassed strtod test...\n");
//
//------------------------- Test getline() ---------------------------
//
//
close( "test_getl" );
x = getline( "test_getl" );
if ( x.[1] != 123.456 ) { error(); }
if ( x.[2] != -123.456 ) { error(); }
if ( x.[3] != 123.456 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != .123 ) { error(); }
if ( x.[2] != -.123 ) { error(); }
if ( x.[3] != .123 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != 123 ) { error(); }
if ( x.[2] != -123 ) { error(); }
if ( x.[3] != 123 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != 1e6 ) { error(); }
if ( x.[2] != -1e6 ) { error(); }
if ( x.[3] != 1e6 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != 1e5 ) { error(); }
if ( x.[2] != -1e5 ) { error(); }
if ( x.[3] != 1e5 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != 123.456e3 ) { error(); }
if ( x.[2] != -123.456e3 ) { error(); }
if ( x.[3] != 123.456e3 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != 123.456e3 ) { error(); }
if ( x.[2] != -123.456e3 ) { error(); }
if ( x.[3] != 123.456e3 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != 123.456e-3 ) { error(); }
if ( x.[2] != -123.456e-3 ) { error(); }
if ( x.[3] != 123.456e-3 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != .123e3 ) { error(); }
if ( x.[2] != -.123e3 ) { error(); }
if ( x.[3] != .123e3 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != 123e3 ) { error(); }
if ( x.[2] != -123e3 ) { error(); }
if ( x.[3] != 123e3 ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != "123abc" ) { error(); }
if ( x.[2] != "abc123" ) { error(); }
if ( x.[3] != "123.abc" ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != "quoted string" ) { error(); }
if ( x.[2] != "q string with escapes \n \t \" " ) { error(); }
x = getline( "test_getl" );
if ( x.[1] != "quoted string" ) { error(); }
if ( x.[2] != 1.23e3 ) { error(); }
if ( x.[3] != 100 ) { error(); }
if ( x.[4] != "q string with escapes \n \t \" " ) { error(); }
if ( x.[5] != 200 ) { error(); }
x = getline ("test_getl", 0);
if (type (x) != "string") { error (); }
// Also check out strsplt while we are here...
if (length (strsplt(x)) != 79) { error (); }
if (length (strsplt(x, 13)) != 6) { error (); }
if (length (strsplt(x,[13,12,14,13,15,11] )) != 6) { error (); }
if (length (strsplt(x, ".")) != 7) { error (); }
printf("\tpassed getline() test...\n");
//
//---------------------- Test readb()/writeb() --------------------
//
a = rand (5,5);
z = rand(3,3) + rand(3,3)*1j;
strm = what()[1:5;1:5];
pi = 4*atan(1.0);
sc = 2*pi;
str = "this is a sample string\ttab";
l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
writeb ("jnk_rb", a, z, strm, sc, str, l);
# Set aside matrices for later tests
check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
clear (a, z, strm, sc, str, l);
close ("jnk_rb");
readb ("jnk_rb");
#
# Now do checks
#
if (!all (all (a == check.a))) { error (); }
if (!all (all (z == check.z))) { error (); }
if (!all (all (strm == check.strm))) { error (); }
if (sc != check.sc) { error (); }
if (str != check.str) { error (); }
if (length (l) != 5) { error (); }
if (!all (all (l.a == check.a))) { error (); }
if (!all (all (l.z == check.z))) { error (); }
if (!all (all (l.strm == check.strm))) { error (); }
if (l.sc != check.sc) { error (); }
if (l.str != check.str) { error (); }
printf("\tpassed binary I/O test...\n");
//
//---------------------- Test read()/write() --------------------
//
a = rand (5,5);
z = rand(3,3) + rand(3,3)*1j;
strm = what()[1:5;1:5];
pi = 4*atan(1.0);
sc = 2*pi;
str = "this is a sample string\ttab";
l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
write ("jnk_ra", a, z, strm, sc, str, l);
# Set aside matrices for later tests
check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
clear (a, z, strm, sc, str, l);
close ("jnk_ra");
read ("jnk_ra");
#
# Now do checks
#
if (a.nr != 5 || a.nc != 5) { error (); }
if (z.nr != 3 || z.nc != 3) { error (); }
if (strm.nr != 5 || strm.nc != 5) { error (); }
if (str != check.str) { error (); }
if (length (l) != 5) { error (); }
if (l.a.nr != 5 || l.a.nc != 5) { error (); }
if (l.z.nr != 3 || l.z.nc != 3) { error (); }
if (l.strm.nr != 5 || l.strm.nc != 5) { error (); }
if (l.str != check.str) { error (); }
printf("\tpassed ascii I/O test...\n");
//
//-------------------------- Test eval () ------------------------------
//
if (1 + 2 != eval("1 + 2")) { error ("eval() error"); }
x = function (s, a) { return eval(s); };
str = "yy = 2 + x(\"2*a\",3)";
z = eval(str);
if (z != yy) { error ("eval() error"); }
printf("\tpassed eval test...\n");
//-------------------------- ------------------------------
printf ("\ttest builtin function correctness...\n");
//-------------------------- ------------------------------
//
//-------------------------- Test abs () ------------------------------
//
A = rand(5,5);
T = ( A == abs (A));
if (!all (all (A))) { error ("abs() incorrect"); }
printf("\tabs()");
//
//-------------------------- Test max () ------------------------------
//
A = [1,10,100;2,20,200;1,2,3];
B = A./2;
ZA = A + rand (3,3)*A*1i;
ZB = B + rand (3,3)*B*1i;
if (!all (max (A) == [2,20,200])) { error( "max() incorrect"); }
if (max (max(A)) != 200) { error ("max() incorrect"); }
if (any (any (max (A, B) != A))) { error (); }
if (any (any (max (B, A) != max (A, B)))) { error (); }
if (any (any (max (ZB, ZA) != max (ZA, ZB)))) { error (); }
if (any (any (max (B, ZA) != max (ZA, B)))) { error (); }
if (any (any (max (ZB, A) != max (A, ZB)))) { error (); }
printf("\tmax()");
//
//-------------------------- Test min () ------------------------------
//
if (!all (min (A) == [1,2,3])) { error( "min() incorrect"); }
if (min (min(A)) != 1) { error ("min() incorrect"); }
if (any (any (min (A, B) != B))) { error (); }
if (any (any (min (B, A) != min (A, B)))) { error (); }
if (any (any (min (ZB, ZA) != min (ZA, ZB)))) { error (); }
if (any (any (min (B, ZA) != min (ZA, B)))) { error (); }
if (any (any (min (ZB, A) != min (A, ZB)))) { error (); }
printf("\tmin()");
//
//-------------------------- Test maxi () -----------------------------
//
if (!all (maxi (A) == [2,2,2])) { error( "maxi() incorrect"); }
if (maxi (maxi(A)) != 1) { error ("maxi() incorrect"); }
printf("\tmaxi()");
//
//-------------------------- Test mini () -----------------------------
//
if (!all (mini (A) == [1,3,3])) { error( "mini() incorrect"); }
if (mini (mini(A)) != 1) { error ("mini() incorrect"); }
printf("\tmini()");
//
//-------------------------- Test floor () ----------------------------
//
if (floor (1.9999) != 1) { error ("floor() output incorrect"); }
if (!all (all (floor ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
error ("floor output incorrect");
}
printf("\tfloor()");
//
//-------------------------- Test ceil () 0----------------------------
//
if (ceil (1.9999) != 2) { error ("ceil() output incorrect"); }
if (!all (all (ceil ([1.99,1.99;2.99,2.99]) == [2,2;3,3]))) {
error ("ceil output incorrect");
}
printf("\tceil()\n");
//
//-------------------------- Test round () ----------------------------
//
if (round (1.8) != 2) { error ("round() output incorrect"); }
if (round (1.4) != 1) { error ("round() output incorrect"); }
if (!all (all (round ([1.99,1.99;2.4,2.4]) == [2,2;2,2]))) {
error ("round output incorrect");
}
printf("\tround()");
//
//-------------------------- Test sum () ------------------------------
//
S = [1:4; 4:7; 8:11];
if (sum (S[1;]) != 10) { error ("sum() incorrect"); }
if (sum (S[3;]) != 38) { error ("sum() incorrect"); }
if (!all (all (sum (S) == [13,16,19,22]))) { error ("sum() incorrect"); }
printf("\tsum()");
//
//-------------------------- Test cumsum () ------------------------------
//
S = [1:4; 4:7; 8:11];
if (any(cumsum (S[1;]) != [1,3,6,10])) { error ("cumsum() incorrect"); }
if (any(cumsum (S[;3]) != [3;9;19])) { error ("cumsum() incorrect"); }
if (!all (all (cumsum (S) == [1,2,3,4;5,7,9,11;13,16,19,22])))
{
error ("cumsum() incorrect");
}
printf("\tcumsum()");
//
//-------------------------- Test prod () ------------------------------
//
S = [1:4; 4:7; 8:11];
if (prod (S[1;]) != 24) { error ("prod() incorrect"); }
if (prod (S[;3]) != 180) { error ("prod() incorrect"); }
if (!all (all (prod (S) == [32,90,180,308])))
{
error ("prod() incorrect");
}
printf("\tprod()");
//
//-------------------------- Test cumprod () ------------------------------
//
S = [1:4; 4:7; 8:11];
if (any(cumprod (S[1;]) != [1,2,6,24])) { error ("cumprod() incorrect"); }
if (any(cumprod (S[;3]) != [3;18;180])) { error ("cumprod() incorrect"); }
if (!all (all (cumprod (S) == [1,2,3,4;4,10,18,28;32,90,180,308])))
{
error ("cumprod() incorrect");
}
printf("\tcumprod()\n");
//
//-------------------------- Test int () ------------------------------
//
if (int (1.9999) != 1) { error ("int() output incorrect"); }
if (!all (all (int ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
error ("int() output incorrect");
}
printf("\tint()");
//
//-------------------------- Test mod () ------------------------------
//
if (mod (1,1) != 0) { error ("mod() output incorrect"); }
if (mod (4,2) != 0) { error ("mod() output incorrect"); }
if (mod (3,2) != 1) { error ("mod() output incorrect"); }
if (mod (5,3) != 2) { error ("mod() output incorrect"); }
printf("\tmod()");
//
//-------------------------- Test find () ------------------------------
//
if (find ([0,1]) != 2) { error ("find() output incorrect"); }
if (find ([1,0]) != 1) { error ("find() output incorrect"); }
if (find ([0,1+1i]) != 2) { error ("find() output incorrect"); }
if (find ([1+0i,0]) != 1) { error ("find() output incorrect"); }
if (find ([0+1i,0]) != 1) { error ("find() output incorrect"); }
printf("\tfind()\n");
//-------------------------- ------------------------------
printf ("\ttest builtin function operation...\n");
//-------------------------- ------------------------------
//
// At least use most of the builtins....
//
A = rand(5,5);
Z = rand(5,5) + rand(5,5)*1j;
S = what()[3:6;2:4];
abs (A);
abs ([A]);
abs (3.14);
abs (Z);
abs ([Z]);
abs (3.14j);
printf ("\tabs");
acos (A);
acos ([A]);
acos (3.14);
acos (Z);
acos ([Z]);
acos (3.14j);
printf ("\tacos");
all (A);
all ([A]);
all (3.14);
all (Z);
all ([Z]);
all (3.14j);
printf ("\tall");
any (A);
any ([A]);
any (3.14);
any (Z);
any ([Z]);
any (3.14j);
printf ("\tany");
asin (A);
asin ([A]);
asin (3.14);
asin (Z);
asin ([Z]);
asin (3.14j);
printf ("\tasin");
atan (A);
atan ([A]);
atan (3.14);
atan (Z);
atan ([Z]);
atan (3.14j);
printf ("\tatan");
atan2 (A,A);
atan2 ([A],[A]);
atan2 (3.14,3.14);
#atan2 (Z,Z);
#atan2 ([Z],[Z]);
#atan2 (3.14j,3.14j);
printf ("\tatan2");
printf ("\n");
cd ("<RLaB$dir>");
printf ("\tcd");
ceil (A);
ceil ([A]);
ceil (3.14);
ceil (Z);
ceil ([Z]);
ceil (3.14j);
printf ("\tceil");
class (A);
class ([A]);
class (3.14);
class ("string");
class (S);
class ([S]);
class (<< rand(3,3); rand(4,4) >>);
class (Z);
class ([Z]);
class (3.14j);
printf ("\tclass");
clear (rand(3,3));
clear ([rand(3,3)]);
printf ("\tclear");
conj (A);
conj ([A]);
conj (3.14);
conj (Z);
conj ([Z]);
conj (3.14j);
printf ("\tconj");
cos (A);
cos ([A]);
cos (3.14);
cos (Z);
cos ([Z]);
cos (3.14j);
printf ("\tcos");
diag (A);
diag ([A]);
diag (3.14);
diag (Z);
diag ([Z]);
diag (3.14j);
printf ("\tdiag");
printf ("\n");
exist (A);
exist (Z);
printf ("\texist");
exp (A);
exp ([A]);
exp (3.14);
exp (Z);
exp ([Z]);
exp (3.14j);
printf ("\texp");
find (A);
find ([A]);
find (3.14);
find (Z);
find ([Z]);
find (3.14j);
printf ("\tfind");
floor (A);
floor ([A]);
floor (3.14);
floor (Z);
floor ([Z]);
floor (3.14j);
printf ("\tfloor");
format (1,1);
format ([2],[3]);
format ();
format ([3], 3);
format ();
printf ("\tformat");
imag (A);
imag ([A]);
imag (3.14);
imag (Z);
imag ([Z]);
imag (3.14j);
printf ("\timag");
inf ();
printf ("\tinf");
printf ("\n");
int (A);
int ([A]);
int (3.14);
int (Z);
int ([Z]);
int (3.14j);
printf ("\tint");
issymm (A);
issymm ([A]);
issymm (3.14);
issymm (Z);
issymm ([Z]);
issymm (3.14j);
printf ("\tissymm");
length (A);
length ([A]);
length (3.14);
length (Z);
length ([Z]);
length (3.14j);
length (3.14);
length (S);
length (<< rand(2,2); rand(10,10); 3.14 >>);
printf ("\tlength");
log (A);
log ([A]);
log (3.14);
log (Z);
log ([Z]);
log (3.14j);
printf ("\tlog");
log10 (A);
log10 ([A]);
log10 (3.14);
log10 (Z);
log10 ([Z]);
log10 (3.14j);
printf ("\tlog10");
max (A);
max ([A]);
max (3.14);
max (Z);
max ([Z]);
max (3.14j);
printf ("\tmax");
maxi (A);
maxi ([A]);
maxi (3.14);
maxi (Z);
maxi ([Z]);
maxi (3.14j);
printf ("\tmaxi");
printf ("\n");
members (<< rand(2,2); rand(); rand(3,3) >>);
printf ("\tmembers");
min (A);
min ([A]);
min (3.14);
min (Z);
min ([Z]);
min (3.14j);
printf ("\tmin");
mini (A);
mini ([A]);
mini (3.14);
mini (Z);
mini ([Z]);
mini (3.14j);
printf ("\tmini");
mod (A,A);
mod ([A],A);
mod (3.14,2);
mod (Z,Z);
mod ([Z],Z);
mod (3.14j,2);
printf ("\tmod");
nan ();
ones(3,3);
ones([4,4]);
printf ("\tones");
prod (A);
prod ([A]);
prod (3.14);
prod (Z);
prod ([Z]);
prod (3.14j);
printf ("\tprod");
real (A);
real ([A]);
real (3.14);
real (Z);
real ([Z]);
real (3.14j);
printf ("\treal");
printf ("\n");
reshape (A, A.nr*A.nc, 1);
reshape (A, [A.nr*A.nc], [1]);
reshape (Z, Z.nr*Z.nc, 1);
reshape (Z, [Z.nr*Z.nc], [1]);
printf ("\treshape");
round (A);
round ([A]);
round (3.14);
round (Z);
round ([Z]);
round (3.14j);
printf ("\tround");
show(A);
show(Z);
show([A]);
show([Z]);
printf ("\tshow");
sign (A);
sign ([A]);
sign (3.14);
sign (Z);
sign ([Z]);
sign (3.14j);
printf ("\tsign");
sin (A);
sin ([A]);
sin (3.14);
sin (Z);
sin ([Z]);
sin (3.14j);
printf ("\tsin");
size (A);
size ([A]);
size (3.14);
size (Z);
size ([Z]);
size (3.14j);
printf ("\tsize");
sizeof (A);
sizeof ([A]);
sizeof (3.14);
sizeof (Z);
sizeof ([Z]);
sizeof (3.14j);
printf ("\tsizeof");
printf ("\n");
sort (A);
sort ([A]);
sort (3.14);
sort (Z);
sort ([Z]);
sort (3.14j);
printf ("\tsort");
sqrt (A);
sqrt ([A]);
sqrt (3.14);
sqrt (Z);
sqrt ([Z]);
sqrt (3.14j);
printf ("\tsqrt");
srand ();
srand ("clock");
srand (SEED);
printf ("\tsrand");
strsplt (S);
printf ("\tstrsplt");
strtod ("string");
strtod (S);
strtod ("1.2");
strtod (["1.2", "1e3"]);
printf ("\tstrtod");
sum (A);
sum ([A]);
sum (3.14);
sum (Z);
sum ([Z]);
sum (3.14j);
printf ("\tsum");
tan (A);
tan ([A]);
tan (3.14);
tan (Z);
tan ([Z]);
tan (3.14j);
printf ("\ttan");
printf ("\n");
tmpnam ();
printf ("\ttmpnam");
type (A);
type ([A]);
type (Z);
type ([Z]);
type (3.14);
type (S);
type ("string");
printf ("\ttype");
zeros (3,3);
zeros ([3], [3]);
zeros ([3,3]);
printf ("\tzeros");
printf ("\n");
srand(SEED);
rand("default");
//
//-------------------------- Test rand () -----------------------------
//
rand ("normal", 5, 1);
xrand = rand(4000,1);
mean = function(x)
{
local(m);
m = size (x)[1];
if( m == 1 )
{
m = size (x)[2];
}
return sum( x ) / m;
};
std = function(x)
{
local(i, m, s);
if(class(x) != "num") { error("std() requires NUMERICAL input"); }
m = x.nr;
if( m == 1 )
{
return sqrt( sum( (x - mean(x)) .^ 2 ) / (x.nc - 1) );
else
for( i in 1:x.nc) {
s[i] = sqrt( sum( (x[;i] - mean(x[;i])) .^ 2 ) / (x.nr - 1) );
}
return s;
}
};
if (!(mean (xrand) > 4.9 && mean (xrand) < 5.1))
{ error ("error in random"); }
if (!(std (xrand) > 0.9 && std (xrand) < 1.1))
{ error ("error in random"); }
printf("\tpassed rand test...\n");
rand("default");
//
//-------------------------- Test norm () -----------------------------
//
tn = [1,2,3,4;2,1,2,3;3,2,1,2;4,3,2,1 ];
if (norm(tn,"m") != 4 ) { error ("incorrect norm computation"); }
if (norm(tn,"1") != 10) { error ("incorrect norm computation"); }
if (norm(tn,"i") != 10) { error ("incorrect norm computation"); }
printf("\tpassed norm test...\n");
//
//-------------------------- Test qr () ------------------------------
//
a = ohess(4);
qa = qr (a);
if (max (max (abs (qa.q*qa.r - a)))/(X*norm (a)*a.nr) > eps)
{ error ("possible qr() problems"); }
z = ohess (4) + ohess(4)*1i;
qz = qr (z);
if (max (max (abs (qz.q*qz.r - z)))/(X*norm (z)*z.nr) > eps)
{ error ("possible qr() problems"); }
printf("\tpassed qr test...\n");
//
//-------------------------- Test schur () ----------------------------
//
a = randsvd (10, 10);
sa = schur (a);
if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
{ error ("possible schur() problems"); }
z = rand (4,4) + rand(4,4)*1i;
sz = schur (z);
if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
{ error ("possible schur() problems"); }
a = randsvd (10, -10);
sa = schur (a);
if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
{ error ("possible schur() problems"); }
z = rand (4,4) + rand(4,4)*1i;
sz = schur (z);
if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
{ error ("possible schur() problems"); }
printf("\tpassed schur test...\n");
//
//-------------------------- Test schord () ----------------------------
//
s2a = schord (sa, 2, 4);
if (max (max (abs (s2a.z*s2a.t*s2a.z' - a)))/(X*norm (a)*a.nr) > eps)
{
error ("possible schord() problems");
}
s2z = schord (sz, 3, 1);
if (max (max (abs (s2z.z*s2z.t*s2z.z' - z)))/(X*norm (z)*z.nr) > eps)
{
error ("possible schord() problems");
}
printf("\tpassed schord test...\n");
//
//-------------------------- Test chol () -----------------------------
//
c = lehmer(10);
u = chol (c);
if (max (max (abs (u'*u - c)))/(X*norm (c)*c.nr) > eps)
{
error ("possible chol() problems");
}
cz = lehmer(10) + lehmer(10)*1j;
cz = symm(cz);
uz = chol (cz);
if (max (max (abs (uz'*uz - cz)))/(X*norm (cz)*cz.nr) > eps)
{
error ("possible chol() problems");
}
printf("\tpassed chol test...\n");
//
//--------------------------- Test inv () --------------------------------
//
a = randsvd(10,10);
b = ones(10,1);
x = inv(a) * b;
if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
{
printf ("\tThe condition # of a: %d\n", 1/rcond (a));
printf ("\tA*X - B:\n");
abs (a*x - b)
error ("possible inv() problems\n");
}
az = randsvd(10,10) + randsvd(10,10)*1j;
bz = rand(10,1) + rand(10,1)*1j;
xz = inv (az)*bz;
if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
{
printf ("\tThe condition # of z: %d\n", 1/rcond (az));
printf ("\tA*X - B:\n");
abs (az*xz - bz)
error ("possible inv() problems\n");
}
printf("\tpassed inv test...\n");
//
//-------------------------- Test solve () -----------------------------
//
//
// Real - General case
//
a = randsvd(10,10);
b = ones(10,1);
x = solve (a,b);
if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
{
printf ("\tThe condition # of a: %d\n", 1/rcond (a));
printf ("\tA*X - B:\n");
abs (a*x - b)
error ("possible solve() problems\n");
}
//
// Real - Symmetric case
//
s = symm (randsvd(10,10));
b = ones(10,1);
x = solve (s,b);
if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
{
printf ("\tThe condition # of s: %d\n", 1/rcond (s));
printf ("\tA*X - B:\n");
abs (s*x - b)
error ("possible solve() problems\n");
}
//
// Complex - General case
//
az = randsvd(10,10) + randsvd(10,10)*1j;
bz = rand(10,1) + rand(10,1)*1j;
xz = solve (az,bz);
if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
{
printf ("\tThe condition # of z: %d\n", 1/rcond (az));
printf ("\tA*X - B:\n");
abs (az*xz - bz)
error ("possible solve() problems\n");
}
//
// Complex - Symmetric case
//
sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
bz = rand(10,1) + rand(10,1)*1j;
xz = solve (sz,bz);
if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
{
printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
printf ("\tA*X - B:\n");
abs (sz*xz - bz)
error ("possible solve() problems\n");
}
printf("\tpassed solve test...\n");
//
//-------------------------- Test factor() / backsub() -----------------------------
//
//
// Real - General case
//
a = randsvd(10,10);
b = ones(10,1);
f = factor (a);
x = backsub (f,b);
if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
{
printf ("\tThe condition # of a: %d\n", 1/rcond (a));
printf ("\tA*X - B:\n");
abs (a*x - b)
error ("possible factor/backsub problems\n");
}
//
// Real - Symmetric case
//
s = symm (randsvd(10,10));
f = factor (s);
x = backsub (f,b);
if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
{
printf ("\tThe condition # of s: %d\n", 1/rcond (s));
printf ("\tA*X - B:\n");
abs (s*x - b)
error ("possible factor/backsub problems\n");
}
s = symm (randsvd(10,10));
f = factor (s, "s");
x = backsub (f,b);
if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
{
printf ("\tThe condition # of s: %d\n", 1/rcond (s));
printf ("\tA*X - B:\n");
abs (s*x - b)
error ("possible factor/backsub problems\n");
}
//
// Complex - General case
//
az = randsvd(10,10) + randsvd(10,10)*1j;
bz = rand(10,1) + rand(10,1)*1j;
fz = factor (az);
xz = backsub (fz,bz);
if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
{
printf ("\tThe condition # of z: %d\n", 1/rcond (az));
printf ("\tA*X - B:\n");
abs (az*xz - bz)
error ("possible factor/backsub problems\n");
}
az = randsvd(10,10) + randsvd(10,10)*1j;
bz = rand(10,1) + rand(10,1)*1j;
fz = factor (az, "g");
xz = backsub (fz,bz);
if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
{
printf ("\tThe condition # of z: %d\n", 1/rcond (az));
printf ("\tA*X - B:\n");
abs (az*xz - bz)
error ("possible factor/backsub problems\n");
}
//
// Complex - Symmetric case
//
sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
fz = factor(sz);
xz = backsub (fz,bz);
if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
{
printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
printf ("\tA*X - B:\n");
abs (sz*xz - bz)
error ("possible factor/backsub problems\n");
}
sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
fz = factor(sz, "s");
xz = backsub (fz,bz);
if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
{
printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
printf ("\tA*X - B:\n");
abs (sz*xz - bz)
error ("possible factor/backsub problems\n");
}
printf("\tpassed factor/backsub test...\n");
//
//------------------------------ Test lu() ---------------------------------
//
static (swap);
lu = function ( A )
{
local (i, l, u, pvt, x)
if (A.nr != A.nc) { error ("lu() requires square A"); }
x = factor (A, "g"); // Do the factorization
//
// Now create l, u, and pvt from lu and pvt.
//
l = tril (x.lu, -1) + eye (size (x.lu));
u = triu (x.lu);
pvt = eye (size (x.lu));
//
// Now re-arange the columns of pvt
//
for (i in 1:max (size (x.lu)))
{
pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
}
return << l = l; u = u; pvt = pvt >>;
};
//
// In vector V, swap elements I, J
//
swap = function ( V, I, J )
{
local (v, tmp);
v = V;
tmp = v[I];
v[I] = v[J];
v[J] = tmp;
return v;
};
a = randsvd(10,10);
lua = lu (a);
if (max (max (abs(a - lua.pvt*lua.l*lua.u)))/(X*norm (a)*a.nr) > eps)
{
printf ("\tThe condition # of a: %d\n", 1/rcond (a));
printf ("\tA - p*l*u:\n");
abs (a - lua.pvt*lua.l*lua.u)
error ("possible lu()/factor() problems\n");
}
//
// Real
az = randsvd(10,10) + randsvd(10,10)*1j;
luz = lu (az);
if (max (max (abs (az - luz.pvt*luz.l*luz.u)))/(X*norm (az)*az.nr) > eps)
{
printf ("\tThe condition # of z: %d\n", 1/rcond (az));
printf ("\tA - p*l*u:\n");
abs (az - luz.ovt*luz.l*luz.u)
error ("possible lu()/factor()() problems\n");
}
printf("\tpassed lu/factor test...\n");
//
//-------------------------- Test svd () -----------------------------
//
a = randsvd(10,10);
s = svd (a);
if (max (max (abs (s.u*diag(s.sigma)*s.vt - a)))/(X*norm (a)*a.nr) > eps)
{
error ("possible svd() problems");
}
z = randsvd(10,10) + rand(10,10)*1j;
sz = svd (z);
if (max (max (abs (sz.u*diag(sz.sigma)*sz.vt - z)))/(X*norm (z)*z.nr) > eps)
{
error ("possible svd() problems");
}
printf("\tpassed svd test...\n");
//
//-------------------------- Test hess () -----------------------------
//
a = randsvd(10,10);
h = hess (a);
if (max (max (abs (h.p*h.h*h.p' - a)))/(X*norm (a)*a.nr) > eps)
{
error ("possible hess() problems");
}
z = randsvd(10,10) + randsvd(10,10)*1j;
hz = hess (z);
if (max (max (abs (hz.p*hz.h*hz.p' - z)))/(X*norm (z)*z.nr) > eps)
{
error ("possible hess() problems");
}
printf("\tpassed hess test...\n");
//
//-------------------------- Test lyap () ------------------------------
//
lyap = function ( A, B, C )
{
local (A, B, C)
if (!exist (B))
{
B = A'; // Solve the special form: A*X + X*A' = -C
}
if ((A.nr != A.nc) || (B.nr != B.nc) || (C.nr != A.nr) || (C.nc != B.nr)) {
error ("Dimensions do not agree.");
}
//
// Schur decomposition on A and B
//
sa = schur (A);
sb = schur (B);
//
// transform C
//
tc = sa.z' * C * sb.z;
X = sylv (sa.t, sb.t, tc);
//
// Undo the transformation
//
X = sa.z * X * sb.z';
return X;
};
a = randsvd (10,10);
b = rand (10,10);
c = rand (10,10);
x = lyap (a, b, c);
if (max (max (abs (a*x + x*b + c)))/(X*norm(c)*norm(a)*norm(b)) > eps)
{
error ("possible problems with lyap() or sylv()");
}
printf("\tpassed lyap test...\n");
//
//-------------------------- Test eig () ------------------------------
//
trace = function(m)
{
local(i, tr);
if(m.class != "num") {
error("must provide NUMERICAL input to trace()");
}
tr = 0;
for(i in 1:min( [m.nr, m.nc] )) {
tr = tr + m[i;i];
}
return tr;
};
eye = function( m , n )
{
local(i, N, new);
if (!exist (n))
{
if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
new = zeros (m[1], m[2]);
N = min ([m[1], m[2]]);
else
if (class (m) == "string" || class (n) == "string") {
error ("eye(), string arguments not allowed");
}
if (max (size (m)) == 1 && max (size (n)) == 1)
{
new = zeros (m[1], n[1]);
N = min ([m[1], n[1]]);
else
error ("matrix arguments to eye() must be 1x1");
}
}
for(i in 1:N)
{
new[i;i] = 1.0;
}
return new;
};
//
// Standard eigenvalue problem
//
a = randsvd(10,10);
ta = trace (a);
sa = symm (a);
tsa = trace (sa);
z = randsvd(10,10) + randsvd(10,10)*1i;
tz = trace (z);
sz = symm (z);
tsz = trace (sz);
tol = 1.e-6;
if (!(ta < sum(eig(a).val) + tol && ta > sum(eig(a).val) - tol))
{
error ("error in eig");
}
if (!(tsa < sum(eig(sa).val) + tol && tsa > sum(eig(sa).val) - tol))
{
error ("error in eig");
}
if (abs(tz)+tol < abs(sum(eig(z).val)) && abs(tz)+tol > abs(sum(eig(z).val)))
{
error ("error in eig");
}
if (abs(tsz)+tol < abs(sum(eig(sz).val)) && abs(tsz) > abs(sum(eig(sz).val)))
{
error ("error in eig");
}
//
// Generalized eigenvalue problem
//
b = randsvd(10,10);
sb = symm (b) + eye(size(b))*3;
tb = trace (b);
tsb = trace (sb);
zb = randsvd(10,10) + randsvd(10,10)*1i;
szb = symm (zb) + eye(size(zb))*3;
tzb = trace (zb);
tszb = trace (szb);
eig(a,b); // not sure of a good way to check these yet
eigs(sa,sb);
eigs(sa,sb);
eig(z, zb);
eigs(sz, szb);
eigs(sz, szb);
printf("\tpassed eig test...\n");
//
//-------------------------- Test fft () -----------------------------
//
if (100 != fft(ones(100,1))[1]) { error ("error in fft()"); }
printf("\tpassed fft test...\n");
//
//------------------------- Fibonacci Test -------------------------
//
// Calculate Fibonacci numbers
//
i=1;
while ( i < 2 ) {
i=i+1;
a=0; b=1;
while ( b < 10000 ) {
c = b;
b = a+b;
a = c;
}
}
if ( b != 10946 ) {
error("failed fibonacci test");
else
printf("\tpassed fibonacci test...\n");
}
//
//------------------------- Factorial Test -------------------------
//
fac = function(a)
{
if(a <= 1)
{
return 1
else
return a*$self(a-1)
}
};
if(fac(10) != 3628800)
{
error(); else printf("\tpassed factorial test...\n");
}
//
//--------------------------- ACK Test ----------------------------
//
ack = function(a, b)
{
if(a == 0) { return b + 1; }
if(b == 0) { return $self(a-1, 1); }
return $self(a-1, $self(a, b-1));
};
if(ack(2,2) != 7)
{
error("error in ack() test");
else
printf("\tpassed ACK test...\n");
}
//
//------------------------- Prime Test -----------------------------
//
// An example that finds all primes less than limit
//
primes = function (limit)
{
local(prime, cnt, i, j, k);
i = 1; cnt = 0;
for(k in 2:limit)
{
j = 2;
while(mod(k,j) != 0)
{
j++;
}
if(j == k) // Found prime
{
cnt++;
prime[i;1] = k;
i++;
}
}
return prime;
};
if(max(size(primes(100))) != 25)
{
error("error in prime test");
else
printf("\tpassed prime test...\n");
}
//
//--------------------------- Fib Min Test -----------------------------
//
// fibmin() will minimize an arbitrary function
// in 1D using Fibonacci search
f065 = function(x)
{
return (x - 0.65) * (x - 0.65);
};
fib = function(x)
{
local (n, a, b);
a = 1;
b = 1;
if (x >= 2)
{
n = x - 1;
for (n in n:1:-1)
{
c = b;
b = a + b;
a = c;
n = n - 1;
}
}
return b;
};
// Minimize a 1D function using Fibonacci search
// f = function to minimize
// xlo = lower bound
// xhi = upper bound
// n = number of iterations (the bigger the more accurate)
fibmin = function(f, xlo, xhi, n)
{
local(a, b, x, y, ex, ey, k, lo, hi);
lo = xlo;
hi = xhi;
k = n;
for (k in k:2:-1)
{
a = fib(k - 2) / fib(k);
b = fib(k - 1) / fib(k);
x = lo + (hi - lo) * a;
y = lo + (hi - lo) * b;
ex = f(x);
ey = f(y);
if (ex >= ey)
{
lo = x;
else
hi = y;
}
// printf("%d: (%g %g) %g %g %g %g\n", k, a, b, lo, hi, ex, ey);
}
return (lo + hi) / 2;
};
//
// Simple example using above function to mimize f065. Answer is 0.65
//
x = fibmin(f065, 0, 1, 30); // printf("f(%g)=%g\n", x, f065(x));
if (abs(x - 0.65) > 1e-6)
{
printf("x = %f\n", x);
error("failed fibmin test");
}
printf("\tpassed fibmin test...\n");
//
//--------------------- Nasty Function Test ------------------------
//
printf("\tStarting Nasty Function Test...");
printf("\tthis will take awhile\n");
check = function( a, b, c, d, e, f, g, h )
{
if ( a+b+c+d == e+f+g+h && ...
a^2+b^2+c^2+d^2 == e^2+f^2+g^2+h^2 && ...
a^3+b^3+c^3+d^3 == e^3+f^3+g^3+h^3 )
{
return 1;
else
return 0;
}
};
for(a in 8:10) {
for(b in 7:(a-1)) {
for(c in 6:(b-1)) {
for(d in 5:(c-1)) {
for(e in 4:(d-1)) {
for(f in 3:(e-1)) {
for(g in 2:(f-1)) {
for(h in 1:(g-1)) {
if(check( a, b, c, d, e, f, g, h ) || ...
check( a, e, c, d, b, f, g, h ) || ...
check( a, f, c, d, e, b, g, h ) || ...
check( a, g, c, d, e, f, b, h ) || ...
check( a, h, c, d, e, f, g, b ) || ...
check( a, b, e, d, c, f, g, h ) || ...
check( a, b, f, d, e, c, g, h ) || ...
check( a, b, g, d, e, f, c, h ) || ...
check( a, b, h, d, e, f, g, c ) || ...
check( a, b, c, e, d, f, g, h ) || ...
check( a, b, c, f, e, d, g, h ) || ...
check( a, b, c, g, e, f, d, h ) || ...
check( a, b, c, h, e, f, g, d ) || ...
check( a, e, f, d, b, c, g, h ) || ...
check( a, e, g, d, b, f, c, h ) || ...
check( a, e, h, d, b, f, g, c ) || ...
check( a, f, g, d, e, b, c, h ) || ...
check( a, f, h, d, e, b, g, c ) || ...
check( a, g, h, d, e, f, b, c ) || ...
check( a, b, e, f, c, d, g, h ) || ...
check( a, b, e, g, c, f, d, h ) || ...
check( a, b, e, h, c, f, g, d ) || ...
check( a, b, f, g, e, c, d, h ) || ...
check( a, b, f, h, e, c, g, d ) || ...
check( a, b, g, h, e, f, c, d ) || ...
check( a, e, f, g, e, f, g, h ) || ...
check( a, e, f, h, e, f, g, h ) || ...
check( a, e, g, h, e, f, g, h ) || ...
check( a, f, g, h, e, f, g, h ) ) { cnt++; }
}
}
}
}
}
}
}
}
if(1) { // figure out the value of cnt, and check!
printf("\tpassed nasty function test...\n");
else
error();
}
//
//------------------ Test More Advanced Functions --------------------
//
printf( "\tStarting the lqr/ode test..." );
printf( "\tthis will take awhile\n" );
lqr = function( a, b, q, r )
{
local( k, s,...
m, n, mb, nb, mq, nq,...
e, v, d );
m = size(a)[1]; n = size(a)[2];
mb = size(b)[1]; nb = size(b)[2];
mq = size(q)[1]; nq = size(q)[2];
if ( m != mq || n != nq )
{
fprintf( "stderr", "A and Q must be the same size.\n" );
quit
}
mr = size(r)[1]; nr = size(r)[2];
if ( mr != nr || nb != mr )
{
fprintf( "stderr", "B and R must be consistent.\n" );
quit
}
nn = zeros( m, nb );
// Start eigenvector decomposition by finding eigenvectors of Hamiltonian:
e = eig( [ a, solve(r',b')'*b'; q, -a' ] );
v = e.vec; d = e.val;
index = sort( real( d ) ).ind;
d = real( d[ index ] );
if ( !( d[n] < 0 && d[n+1] > 0 ) )
{
fprintf( "stderr", "Can't order eigenvalues.\n" );
quit
}
chi = v[ 1:n; index[1:n] ];
lambda = v[ (n+1):(2*n); index[1:n] ];
s = -real(solve(chi',lambda')');
k = solve( r, nn'+b'*s );
return << k=k; s=s >>;
};
// Now run a little test problem.
k = 1; m = 1; c = .1;
a = [0 ,1 ,0 , 0;
-k/m, -c/m, k/m, c/m;
0, 0, 0, 1;
k/m, c/m, -k/m, -c/m ];
b = [ 0; 1/m; 0; 0 ];
qxx = diag( [0, 0, 100, 0] );
ruu = [1];
K = lqr( a, b, qxx, ruu ).k;
dot = function( t, x )
{
global (a, b, K)
return (a-b*K)*x + b*K*([1,0,1,0]');
};
x = ode ( dot, 0, 15, [0,0,0,0], .02, 1e-5, 1e-5 );
m = maxi( x[;2] );
if ( (abs( x[m;2] - 1.195 ) > 0.001) || ...
any (abs( x[x.nr;2,4] - 1 ) > 0.001) )
{
printf( "\tfailed***\n" );
else
printf( "\tpassed the lqr/ode test...\n" );
}
printf("Elapsed time = %10.3f seconds\n", toc() );
"FINISHED TESTS"